- BUILD_GPU.md, GPU_ACCELERATION_PLAN.md: relocate into the repo (were
  untracked at the workspace root) so they are version-controlled.
- BUILD_GPU.md: record CPU-build validation results -- energy parity vs
  one-at-a-time xtb (bit-for-bit on a diverse set) and the batched-eigensolver
  vs per-system gate (max |dEPS| = 1.95e-14 eV with up to 28 padding rows
  exercised, PASS <= 1e-6 eV).
- benchmark/gpu_batch_bench.sh: portable A (per-process) vs B (--gpu-batch)
  throughput benchmark. Dev box: 96 GFN0 structures, 44.11s -> 18.35s = 2.40x
  (CPU startup-amortization win; eigensolver speedup comes on top once the GPU
  backend is wired).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Split the GFN0 single point so the validated batched eigensolver drives the
production total energy (energy + properties; gradients stay on per-molecule
peeq, Phase 3). The original peeq (gradient path for opt/MD) is untouched.

- peeq_module.f90: new peeq_build_energy (prologue -> H/S, gradient discarded
  via a scratch buffer; energies exact) and peeq_finish_energy (density, shell
  charges, band energy, repulsion+SRB, total energy, gap, Wiberg). Reuse peeq's
  private builders (ccm_build_SH0, drep_grad, dsrb_grad). New TPeeqEnergyCtx
  carries the small state across the build -> diagonalize -> finish boundary.
- batch_driver.F90: run_batched_energy builds each molecule's H/S, bins by basis
  size, diagonalizes whole buckets in one TBatchedEigensolver call, scatters
  eigenpairs into each wavefunction, finishes energies, and checks the totals
  against the per-molecule peeq reference.

Validated on CPU: batched total energies match per-molecule peeq to ~1e-15 Eh
(diverse set, padding up to 28 rows; and 96 water dimers). xtb unit suite
unaffected (xtb:peeq, xtb:repulsion OK).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

The GPU numerical gate from the plan, run on real hardware (RTX 3050, cc 8.6):
cusolverDnDsygvd reproduces LAPACK dsygvd to ~1e-14 eV on real GFN0 H/S, both
bare and with the --gpu-batch bucket padding -- far inside the 1e-6 gate.

- test/gpu/cusolver_gate.f90: standalone nvfortran program; reads dumped real
  GFN0 (H,S), solves via cuSolver (the exact path in batched_eig.F90) and LAPACK,
  compares eigenvalues. Bare + padded-bucket tests.
- test/gpu/run_gpu_gate.sh: dump (CPU build) -> compile (nvfortran -acc
  -gpu=cc86 -cudalib=cusolver,cublas) -> run on GPU.
- batch_driver.F90: dumpCapturedHS (XTB_DUMP_HS) writes captured real
  eigenproblems for the gate.
- BUILD_GPU.md: GPU gate results + toolchain note.

Toolchain finding: nvfortran 26.3 (LLVM-only) ICEs sizing deferred-length
character components in xtb's dependency tree (toml-f etc.), so full xtb is not
yet buildable with nvfortran -- hence the standalone gate. OpenACC GPU offload
itself works on the RTX 3050. Full-xtb-on-nvfortran is a separate effort.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

- run_gpu_batch now prints an in-place carriage-return progress bar
  ([####----] pct done/total name) as it processes each structure.
- The dev-only passes (capture + cuSolver-gate replay via validateBatchedEig,
  and the batched-energy parity check via run_batched_energy) are now gated
  behind the XTB_GPU_VALIDATE env var, so ordinary screening runs ~2-3x less
  work and the progress bar reflects the real workload.

Signed-off-by: Prasanna163 <prasannakulkarni163@gmail.com>

…ation)

Adds genuine GPU diagonalization to the gfortran-built xtb, sidestepping the
nvfortran blocker (nvfortran cannot compile xtb's dependency tree). The
eigensolver is plain CUDA C compiled by nvcc and linked in.

- src/gpu/gpu_eig.cu: gpu_sygvd_batch -- batched generalized symmetric-definite
  eigensolver via cuSolver (cusolverDnDsygvd); H overwritten with eigenvectors,
  W gets eigenvalues. Callable from Fortran via iso_c_binding.
- meson: -Dgpu_shim=true compiles gpu_eig.cu with nvcc (custom_target) and links
  cuSolver/cudart into the normal gfortran build; defines WITH_GPU_SHIM. New
  options gpu_shim / nvhpc_root / gpu_cuda_ver. Build static (FC=gfortran).
- batch_driver.F90: iso_c_binding interface + gpu_use flag; run_batched_energy
  now fills the results table, shows a progress bar, and routes the per-bucket
  diagonalization to the GPU when --gpu is set (XTB_BATCH_CPU forces CPU for
  benchmarking). run_gpu_batch gains a GPU-primary path.
- main.F90: --gpu flag (implies --gpu-batch; GFN0).

Validated on an RTX 3050: `xtb --gfn 0 --gpu` energies match CPU exactly.
Benchmark (300 real molecules, same path, only diag differs): GPU 32.5 s vs
CPU 38.2 s = 1.18x end-to-end -- diag is GPU-fast but a small slice of each
small-molecule single point; 16-thread xtbfolder (parallel CPU) is still
fastest for small molecules. GPU wins on larger systems.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Signed-off-by: Prasanna163 <prasannakulkarni163@gmail.com>

The single-process --gpu path built each molecule's H/S serially, so the
GPU-accelerated diagonalization was bottlenecked by the serial CPU build/EEQ/
dispersion work. Parallelize the per-molecule build loop with OpenMP:

- Each molecule uses a thread-private TEnvironment (lenv), so the shared error
  log is never raced; control flow keys off ctx%fail/orthog, not the global env.
- Per-thread private gbsa; reductions for the built/fail/orthog counters; the
  progress bar update is in a critical section. paramFile is a fixed-length
  read-only shared string. Nested OpenMP stays disabled, so each molecule's
  internal integral OMP runs serially inside its parallel iteration (no
  oversubscription). OMP_NUM_THREADS controls the number of parallel molecules.

Verified on the RTX 3050 (300 real molecules, GFN0):
- 3 repeated GPU runs byte-identical (deterministic; no race observed)
- GPU(OMP=8) energies == CPU(OMP=8) energies for all 300 (parity OK)
- parallel build 3.0x over serial (12.4 -> 4.2 s); GPU 1.72x over CPU same path

Note: xtb is not generally thread-safe across molecules; this is validated for
the H/C/N/O screening set. Energy + properties only (no gradient).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Signed-off-by: Prasanna163 <prasannakulkarni163@gmail.com>

Code-grounded plan for the self-consistent GPU path: lockstep batched SCF with a
per-molecule active/converged mask. Identifies the scc SCF loop (scc_core.f90),
the refactor (scc -> init/step/final), the batched kernels (diag, dmat GEMM, H1
build, GFN2 AES), 6 validated increments, and honest scope (weeks; GPU wins on
large systems / big batches, not tiny molecules).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Signed-off-by: Prasanna163 <prasannakulkarni163@gmail.com>

…cuBLAS)

Routes the self-consistent GFN1/GFN2 path through the GPU when --gpu is set,
reusing the gfortran build + CUDA-C shim (no nvfortran).

- scc_core: split the monolithic scc into behavior-preserving init/step/final
  seams with an explicit TScfBatchState (factorized overlap, Broyden history,
  GFN2 multipole intermediates, damping/convergence, iteration counters, active
  flag). The SCF iteration's generalized diagonalization and density build route
  to cuSolver/cuBLAS when --gpu is set; CPU fact_solve otherwise.
- gpu_runtime.F90: leaf module owning the gpu_use flag + CUDA eigensolver
  interface + a plain gpu_solve wrapper callable from non-preprocessed core code.
- gpu_eig.cu: persistent CUDA context (handles/workspaces reused across SCF and
  opt cycles); cusolverDnDsygvd, cuBLAS density, CUDA GFN1 isotropic H1 build,
  CUDA shell-Mulliken; CPU fallback on any CUDA error.
- --gpu no longer implies --gpu-batch; works with --sp/--grad/--opt, and with
  --gpu-batch drives the same CUDA path per structure in the multi-file driver.
- meson/cmake: WITH_GPU_SHIM, nvcc custom_target with -arch=sm_<gpu_arch>,
  cusolver/cublas/cudart link. New GFN1/2 GPU gate (test/gpu/) registered.

Validated on RTX 3050 (sm_86), CPU vs GPU, OMP=1:
  GFN1: |dE_sp|=0  max|dG|=9.99e-16  |dE_opt|=0  max|dx|=0
  GFN2: |dE_sp|=0  max|dG|=1.70e-16  |dE_opt|=0  max|dx|=0
Single point, gradient, and loose optimization all match CPU (gate passes).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Signed-off-by: Prasanna163 <prasannakulkarni163@gmail.com>

reportResults writes structure,energy_Eh,gap_eV,status to the path in
XTB_BATCH_CSV when set, so the xtbx folder pipeline can aggregate batch
energies into a summary table. Also ignore the cusolver_gate / gpu_bench
build artifacts.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

The per-call GPU offload copied H/S/C/P host<->device every SCF iteration
(~8 n^2 transfers/iter) and re-uploaded the constant overlap ~3x, so the GPU
was slower than the CPU for small molecules. A resident session
(gpu_scf_open/solve/finish/get_vectors/close in src/gpu/gpu_eig.cu, wrapped in
xtb_gpu_runtime) keeps S -- and for GFN1 H0/matlist/ao2sh -- resident for the
whole SCF; H and the density stay resident between solve and finish, so per
iteration only small vectors cross PCIe (shift/focc in; emo/qsh out; P out for
the host electro()). scc_core opens the session in scc_init, routes solve/
finish in scc_step, fetches the converged eigenvectors once (gradient's
energy-weighted density), and closes in scc_final, with a full CPU fallback.

Gate stays bit-exact (GFN1/GFN2 dE_sp=0, grad ~1e-16, opt dx=0). Benchmarks
(RTX 3050, 8 CPU helper threads): 648-atom water cluster GFN2 95.2s -> 10.1s
(9.5x), GFN1 433.8s -> 13.6s (32x). Small molecules still favor the CPU.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Back up the working one-command xtbx front-end (dispatcher + setup + helpers)
into the repo. Folder runs optimize each compound by default into
results/<name>/ with a live progress bar and an energy table; single large
molecules auto-route to the GPU. Paths inside are machine-specific (documented
in tools/win/README.md).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Pre-relaxes with the O(N) GFN-FF force field before the GFN2/GFN1 optimization
so far-from-minimum inputs (raw/built/docked structures, the protein prototype)
need fewer expensive SCF cycles. Opt-in: already-good crystal geometries don't
benefit (GFN-FF and GFN2 minima differ slightly).

Works around a GFN-FF xtbopt.xyz writer bug on this build (it corrupts the
element-symbol column with arbitrary bytes; coordinates are correct) by
rebuilding the geometry from the ORIGINAL symbols + the GFN-FF-optimized
coordinates, recovered by regex so it is robust to the corruption.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Reproducible CPU-vs-GPU benchmarks (taxol single point, large water cluster,
GFN2 --opt time-per-cycle) plus the measured crossover (~350 atoms) and VRAM
ceiling (~173 bytes/n^2 -> n<=~5500 on 6 GB). README captures the results.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Profiling shows the GFN2 analytical gradient is ~42% of a warm optimization
cycle (build_dSDQH0 integral-derivative contraction dominates), so offloading it
is worthwhile. A pre-existing OpenACC port exists (#ifdef XTB_GPU !$acc regions);
this adds a gpu_acc meson option to drive it and verified gfortran NVPTX offload
works on the RTX 3050.

However the !$acc code is nvfortran+managed-memory-targeted and does NOT build
under gfortran: nested vector parallelism (aespot, fixed -> inner loops seq);
derived types listed whole AND by-component in data clauses (gfortran 'mixed
component access'); acc routine gang in PURE procs; module params in copy
clauses; default(present) regions assuming managed memory. Enabling it is a real
port. SHELVED -- gpu_acc is off by default and marked experimental; production
GPU builds use gpu_shim (cuSolver) only. See GPU_PHASE2_BATCHED_SCF.md Inc 3.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

molecule_to_structure passed mol%sym (character(len=80)) to mctc's
new_structure; on this build the long strings corrupted get_identity_symbol's
id mapping, so EVERY geometry optimization (GFN0/1/2 and GFN-FF) wrote garbage
element symbols for atoms 2..N in xtbopt.xyz. Coordinates were correct and
energies (from stdout) unaffected, so it went unnoticed -- but the written
structures were unusable as input ('Cannot map symbol to atomic number').

Build a clean, canonically-sized symbol array (toSymbol) before constructing
the structure. Verified: GFN2/GFN-FF/GFN0 now write correct symbols (0/113 bad
on taxol) and the output reads back normally.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Plan for the GPU gradient (build_dSDQH0 hotspot, ~42% of a warm opt cycle):
C interface, atom-pair parallelization, and incremental milestones each gated
bit-exact before the next. Reference algorithm is the OpenACC build_dSDQH0_gpu.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…back

First milestone of the CUDA-C analytical-gradient kernel (see
GPU_GRADIENT_KERNEL_PLAN.md). Adds the C entry gpu_build_dsdqh0 (stub returning
non-zero), the xtb_gpu_runtime gpu_grad_dsdqh0 wrapper, and routes scf_module's
gradient through it under --gpu with a CPU build_dSDQH0 fallback. The kernel
body (multipole_grad_3d / dtrf2 / olapp / shiftintg integral-derivative
contraction) is filled in over later milestones.

Verified: build-gpushim compiles and the GFN1/GFN2 gate stays bit-exact
(dE_sp=0, max|dG|~1e-16, opt dx=0) via the CPU fallback -- no regression.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…l derivatives)

Full analytical-gradient kernel for GFN1/GFN2: ports get_grad_multiint
(multipole_grad_3d, horizontal_shift, form_product, olapp, shiftintg), dtrf2
(CAO->SAO), h0scal and dshellPoly to CUDA device functions, with one thread per
atom pair contracting the integral derivatives against P/Pew and the AES
potentials (vs/vd/vq) and atomicAdd into g/sigma/dhdcn. scf_module marshals
hData + basis + potentials (with per-array leading dims) and routes the gradient
through it under --gpu, falling back to CPU build_dSDQH0 on any failure.

VALIDATED machine-precision vs CPU: water GFN1 9.99e-16, water/taxol GFN2
1.7e-16/4.0e-15, H2S GFN2 (d-shells) 3.0e-16. Full GFN1/2 gate (SP+grad+opt)
PASSES bit-exact. (A d-shell dtrf2 index transpose was found and fixed.)
Currently correct but overhead-limited per call; a resident gradient session is
the next optimization.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Two bit-exact kernel optimizations: (1) dtrf2 (CAO->SAO) is a no-op for s/p
shells, so skip the 57 build-transform-store passes unless a d-shell is
involved -- big saving for C/N/O/H; (2) zero-init and process only the touched
[naoj][naoi] Cartesian block instead of the full [6][6]. Cuts per-thread local-
memory traffic. Kernel time at 648 atoms: 301ms -> ~80ms (3.8x); GPU gradient is
now ~8x faster than the CPU build_dSDQH0. Added XTB_GRAD_TRACE phase timing
(upload/kernel/download). Gate PASS (GFN1/2 SP+grad+opt bit-exact); d-shell
(H2S) still 3e-16.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Wall-time accumulators around setvsdq / buildIsoAnisotropicH1 / mmompop /
aniso_electro, printed at SCF end when XTB_AES_TRACE is set. Used to scope the
GFN2 AES-on-GPU work: at 648 atoms the four routines cost ~1.2/1.1/1.5/0.7s
(~44% of the SCF), none dominant; mmompop/buildH1 need dpint/qpint resident.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

One thread per atom; sums the gab3/gab5 multipole interactions over j plus the
per-atom CT correction (dipKernel/quadKernel). Standalone (gab uploaded per
call). scc_step routes via gpu_aes_setvsdq under --gpu, CPU setvsdq fallback.
Gate PASS bit-exact (GFN2 dE_sp=0, grad 5e-17).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Standalone AES kernels were overhead-bound (re-uploading the nat^2 gab matrices
per SCF iter made setvsdq slower than CPU). Add a resident AES context
(gpu_aes_open/close): gab3/gab5/dipKernel/quadKernel/xyz/at stay on the device
for the whole SCF; per iteration only q/dipm/qp cross. setvsdq and aniso_electro
(per-atom CT + atom-pair reduction) now run on GPU. scc opens the context in
scc_init (GFN2), routes both routines, closes in scc_final; CPU fallback intact.

Gate PASS bit-exact (GFN2 dE=0). Per-iter cost: setvsdq 1.18->0.46s,
aniso_electro 0.68->0.19s (648-atom SP, summed over iters).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Completes the GFN2 anisotropic electrostatics on the GPU (Stage B), adding the
two remaining routines on top of the Stage A resident context (setvsdq +
aniso_electro):

  * mmompop  -> aes_mmompop_{off,diag,trace}_kernel: multipole populations
    dipm/qp from the density P, using resident dpint/qpint/S/xyz. qpint read in
    xx,yy,zz,xy,xz,yz order (kj=k+l+1), qp written in lin order; trace removed
    per atom (1.5x scale, subtract 0.5*tr).
  * buildIsoAnisotropicH1 -> aes_buildh1_{mat,dip,quad}_kernel: isotropic+
    anisotropic Fockian from resident H0/S/dpint/qpint + per-iter shellShift/
    vs/vd/vq. matlist spans the full packed triangle incl. the diagonal, so the
    mat kernel fully initializes H; dip/quad ADD afterwards on the same stream.

gpu_aes_open now uploads dpint/qpint/S/H0/aoat2/ao2sh/matlist/mdlst/mqlst and
allocates Pm/Hm/shift scratch so both new routines run from the resident context.

All four AES routines are bit-exact vs the CPU (gate: GFN2 dE=0, max|dG|=1.1e-16,
dx=8e-14).

GPU AES is OPT-IN via XTB_GPU_AES=1 (default off). At pocket scale (hundreds of
atoms) the 8-thread CPU AES is 3-10x faster because these routines are memory/
atomic-bound (mmompop atomicAdd contention; buildH1 H round-trip via host), not
compute-bound. The default --gpu path keeps AES on the CPU while the analytical
gradient and diagonalization stay on the GPU (the proven wins). The gate forces
XTB_GPU_AES=1 so the bit-exact GPU AES path stays validated.

Adds xtb-gpu.sh launcher (sets CUDA runtime env + XTBPATH, appends --gpu).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

nvcc on Windows only accepts MSVC as its host compiler, so its output
can't be linked directly into a MinGW/gfortran build. Instead the shim
is prebuilt out-of-band as a DLL (nvcc+MSVC), and a MinGW-compatible
import library is generated from it via dlltool, which xtb_exe now
links against on host_machine.system() == 'windows'. The Linux nvcc
custom_target path is untouched.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>